Large Fluctuations of the Macroscopic Current in Diffusive Systems: 
A Confirmation of the Additivity Principle 
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Most systems, when pushed out of equilibrium, respond by building up currents of locally- 
conserved observables. Understanding how microscopic dynamics determines the averages and fluc- 
tuations of these currents is one of the main open problems in nonequilibrium statistical physics. 
The additivity principle is a theoretical proposal that allows to compute the current distribution in 
many one-dimensional nonequilibrium systems. Using simulations, we confirm this conjecture in a 
simple and general model of energy transport, both in the presence of a temperature gradient and 
in canonical equilibrium. In particular, we show that the current distribution displays a Gaussian 
regime for small current fluctuations, as prescribed by the central limit theorem, and non-Gaussian 
(exponential) tails for large current deviations, obeying in all cases the Gallavotti- Cohen fluctuation 
theorem. In order to facilitate a given current fluctuation, the system adopts a well-defined temper- 
ature profile different from that of the steady state, and in accordance with the additivity hypothesis 
predictions. System statistics during a large current fluctuation is independent of the sign of the 
current, which implies that the optimal profile (as well as higher-order profiles and spatial correla- 
tions) are invariant upon currenst inversion. We also demonstrate that finite-time joint fluctuations 
of the current and the profile are well described by the additivity functional. These results confirm 
the additivity hypothesis as a general and powerful tool to compute current distributions in many 
nonequilibrium systems. 

PACS numbers: 



I. INTRODUCTION 

Understanding the physics of systems out of equilib- 
rium remains challenging to a large extent, even in the 
simplest setting for which one could expect to make 
significant advances, which is that of a nonequilibrium 
steady state (NESS). Even in this simple situation dif- 
ficulties abound mainly because out of equilibrium the 
dynamics plays a dominant role [HQ. For instance, the 
phase space available to a system in a NESS depends cru- 
cially on the dynamics, resulting in a probability measure 
for microscopic configurations which is not known in gen- 
eral for a NESS, as it will inherit this dependence on the 
dynamics Q . This is in contrast to the equilibrium case, 
where the available phase space is uniquely determined 
by the Hamiltonian and the Gibss distribution provides 
the probability measure for microscopic configurations. 
One can ask however questions on the statistics of the 
macroscopic observables characterizing a NESS, as for in- 
stance the current flowing through the system 0-0] • In 
equilibrium, the fluctuations of macroscopic quantities, 
which are a reflection of the hectic microscopic world, 
are strikingly independent of microscopic details, being 
solely determined by thermodynamic quantities as the 
entropy, free energy, etc. A natural way to seek a macro- 
scopic theory of nonequilibrium phenomena is thus to 
investigate the fluctuations of macroscopic currents. Un- 
veiling the relation between microscopic dynamics and 
current fluctuations has proven to be a difficult task |^- 
[l3j . and up to now only few exactly-solvable cases are 
understood. An important step in this direction has 
been the development of the Gallavotti-Cohen fluctua- 



tion theorem [l2|, [H[ , which relates the probability of for- 
ward and backward currents reflecting the time-reversal 
symmetry of microscopic dynamics. However, we still 
lack a general approach based on few simple principles. 
Recently, Bertini, De Sole, Gabrielli, Jona-Lasinio and 
Landim [J] have introduced a Hydrodynamic Fluctua- 
tion Theory (HFT) to study large dynamic fluctuations 
of diffusive systems. This is a very general approach 
which leads to a hard variational problem whose solu- 
tion remains challengingin most cases. Simultaneously, 
Bodineau and Derrida p-Q have conjectured an addi- 
tivity principle for current fluctuations in one dimension 
which can be readily applied to obtain quantitative pre- 
dictions and, together with HFT, seems to open the door 
to a general theory for nonequilibrium systems. 

In this paper we test in depth the validity of the ad- 
ditivity principle in a simple and very general diffusive 
model. In particular, we investigate the fluctuations of 
the energy current in the one-dimensional (ID) Kipnis- 
Marchioro-Pressuti (KMP) model of heat conduction, 
which represents at a coarse-grained level a large class 
of quasi- ID diffusive systems of technological and theo- 
retical interest for which understanding current statistics 
is of central importance. Our results strongly support 
the validity of the additivity principle to describe cur- 
rent fluctuations in one dimension, both in the presence 
of a temperature gradient (NESS) and in canonical equi- 
librium. In particular, we find that the current distribu- 
tion shows both Gaussian and non-Gaussian regimes, and 
obeys the Gallavotti-Cohen symmetry. The system mod- 
ifies its temperature profile to facilitate a given current 
fluctuation, as predicted by the theory, and this profile 
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(as well as any other higher-order profile and spatial cor- 
relation) turns out to be independent of the sign of the 
current. We also explore physics beyond the additivity 
conjecture by studying the fluctuations of the total en- 
ergy in the system, which exhibit the trace left by correc- 
tions to local equilibrium resulting from the presence of 
weak long-range correlations in the NESS. In addition, 
we extend the additivity hypothesis to study the joint 
fluctuations of the current and the profile. 

The paper is structured as follows. In next section 
we describe the additivity principle from a general per- 
spective. Section UTTl introduces the KMP model in one 
dimension. In section IIVI we report the results of our 
simulations, together with a detailed comparison with 
theoretical predictions. Here we also show evidence of 
structure beyond the additivity scenario. Section [V] in- 
vestigates the joint fluctuations of the current and the 
temperature profile, extending the additivity principle 
to understand these finite-time corrections. Finally, we 
present our conclusions in section I VII and a number of 
appendices describe some technical aspects of the discus- 
sion in the main text. Part of the work reported in this 
paper was presented in a shorter Letter (lj| . 



II. THE ADDITIVITY PRINCIPLE 

The additivity principle (to which we will also refer 
here as BD theory) is a conjecture first proposed by T. 
Bodineau and B. Derrida [5J that enables one to calcu- 
late the fluctuations of the current in ID diffusive systems 
in contact with two boundary thermal baths at different 
temperatures, Tl ^ T R . It is a very general conjecture 
of broad applicability, expected to hold for ID systems 
of classical interacting particles, both deterministic or 
stochastic, independently of the details of the interac- 
tions between the particles or the coupling to the thermal 
reservoirs. The only requirement is that the system at 
hand must be diffusive, i.e. Fourier's law must hold. If 
this is the case, the additivity principle predicts the full 
current distribution in terms of its first two cumulants. 
Equivalently, one may use the same formalism to study 
diffusive particle systems coupled to particle reservoirs at 
the boundaries at different chemical potentials, and obey- 
ing Fick's law, or any other open diffusive system char- 
acterized by a single locally-conserved field. However, in 
this paper we stick for simplicity to the energy-diffusion 
version of the problem. Let Pjv(<7, Tl, T R , t) be the prob- 
ability of observing a time-integrated current Qt = qt 
during a long time t in a system of size N. This proba- 
bility typically obeys a large deviation principle (la , [l6| , 



P J v(< ? ,T L ,T fl ,i)~e+^> T - T «) 



(1) 



where J-Ar(q, T/,, Tr) is the current large-deviation func- 
tion (LDF), such that F N {{q),T L ,T R ) = and F N (q ^ 
(q},T L ,T R ) < 0, with (q) = limt-^ooQt/t- This means 
in particular that current fluctuations away from the av- 
erage are exponentially unlikely in time. The additiv- 



ity principle relates this probability with the product of 
probabilities for sustaining the same current in subsys- 
tems of lengths N — n and n, 



PN(q,T L ,T R ,t) 



max [Pjv- 

T 



t (q,T L ,T,t) P n (q,T,T R ,t)}. 

, ( 2 ) 

The maximization over the contact temperature T can 
be rationalized by writing the above probability as an 
integral over T of the product of probabilities for sub- 
systems and noticing that these should obey also a large 
deviation principle akin to eq. ([1]). Hence a saddle-point 
calculation in the long-i limit leads to @. The additiv- 
ity principle can be then rewritten for the large deviation 
function as 



T N (q,T L ,T R ) 



max \J-n- 
T 



M, T L ,T) + T n (q, T, T R )} 



(3) 

We now may adopt a scaling form FnU^TljTb) = 
N~ 1 Q(Nq,TL,T R ) for the current LDF [5-7], and pro- 
ceed by slicing iteratively the ID system of length N 
into smaller and smaller segments. For small enough seg- 
ments the temperature difference across each of them will 
be small, so for small currents q ~ 0(N~ 1 ) each interval 
can be considered to be close to equilibrium and hence 
exhibits locally-Gaussian fluctuations around the average 
current (given by Fourier's law) at the leading order. In 
this way we obtain in the continuum limit the following 
variational form for Q 0-0] 



Q(q) = - min 




n[T q {x)]T' q {x)\ 
2a[T q (x)} 



■dx 



(4) 



where we dropped the dependence on the baths for con- 
venience. Here n(T) is the thermal conductivity char- 
acterizing Fourier's law, (Q t )/t = -k(T) VT, and cr(T) 
measures current fluctuations in equilibrium (Tl = T R ), 
(Qt)/t = u(T)/N. The optimal temperature profile 
T q (x) derived from Q by functional differentiation obeys 



K 2 [T q {x)} 



dT q (x) 
dx 



{l + 2a[T q (x)}K(q 2 )} , (5) 



where K(q 2 ) is a constant which guarantees the correct 
boundary conditions, T q (Q) = Tl and T q (\) — T R . In 
what follows we assume Tl > T R without loss of gen- 
erality. Equations (@| and (JS|) completely determine the 
current distribution, which is in general non-Gaussian 
(except for very small current fluctuations) and obeys 
the Gallavotti-Cohen symmetry, 



G(-q) = G(q)-£q, 
with £ a constant defined by [5J 



(6) 



£ = 2 



a(T) 



dT. 



Moreover, the optimal profile solution of eq. ([5]) is in- 
dependent of the sign of the current, T q (x) = T- q (x), a 
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rather counter-intuitive result which, together with the 
Gallavotti-Cohen relation, reflects the time-reversal sym- 
metry of microscopic dynamics [HI, [l3j] . 

In the simplest case, when K(q 2 ) is large enough for 
the rhs of eq. ([5]) not to vanish -something that happens 
for currents close to the average, the optimal profile T q (x) 
is monotone and we have (Tj, > Tr) 

Using this expression in eq. (@| leads to 



S(q) 



a(T) 



VI + 2K(q 2 )a(T) 



(8) 

and integrating eq. Q above over the whole interval 
x G [0, 1] we obtain an implicit equation for K(q 2 ), 



k(T) 



t r v/1 + 2K(q 2 )a(T) 



.dT. 



(9) 



In many applications it is interesting to work with the 
Legendre transform of the large deviation function, 



J\ q 



or equivalently /x(A) = N~ 1 [Q(q ) + \q ], with q (X) given 
by d q Q(q ) + \ = 0. By noticing that d q Q{q) = Q/q + Kq, 
it then follows for monotone profiles 



M(A) = - 



K{\) 
N 




y/l + 2K(\)(T(T) 



dT\ , (11) 



where K(X) is now obtained from 



sgn[ 9o (A)] 



v/1 + 2K(\)a(T) 



dT. 



(12) 



and sgn(q) = \q\/q is the sign function. The function 
fx(X) can be viewed as the conjugate potential to G{q), 
with A the parameter conjugate to the current q, a re- 
lation equivalent to the free energy being the Legendre 
transform of the internal energy in thermodynamics, with 
the temperature as conjugate parameter to the entropy. 

When the constant K is negative enough for the rhs 
of eq. (JHJ) to vanish at some point, the resulting opti- 
mal profile T q (x) becomes non- monotone. In this case 
it can be shown [5j that the expressions for G(q) and 
K(q 2 ), or their equivalent formulas in A-space, are just 
the analytic continuation of their monotone-case coun- 
terparts. Appendix [A] shows the particular expressions 
for the current LDF and the associated optimal profile, 
both in the monotone and non-monotonous cases, as de- 
rived when applying this general scheme to the particular 
model of interest in this paper, the Kipnis-Marchioro- 
Presutti (KMP) model of heat conduction [l7j . 



Before continuing with the description of this model, 
it is worth noticing that the additivity principle can be 
better understood within the context of Hydrodynamic 
Fluctuation Theory of Bertini et al. Q, which provides 
a variational principle for the most probable (possibly 
time-dependent) profile responsible of a given current 
fluctuation. The probability of observing a particular his- 
tory of the temperature profile T(x, t) and the rescaled 
current j(x, t) during a macroscopic time is, according to 
HFT0,|3, 

P({T(x,t),j(x,t)}) ~exp(-7V2 t [T ! i]) (13) 

where the functional X t can be written as 

HT,j]=fdr [\ x mT) + 4T( X ,T)]r { x,r)f 
Jo Jo 



2a[T(x,T)} 



(14) 

and where the rescaled current field is related to the tem- 
perature profile via the continuity equation d T T(x, r) + 
d x j(x 1 r) — 0. The large deviation function of the inte- 
grated current is then 



P 



N 



exp 



t 

N 



(10) where Q(q) is related to I t [T,j] via 



Q(q) = lim 



with the constraint 



9=7 



— min I* \T, j] 



j(x,T)dr, 



(15) 



(16) 



(17) 



and T(x, r) and j(x, r) coupled via the above continuity 
equation. Solving this time-dependent problem to obtain 
explicit predictions for the current LDF remains a chal- 
lenge in most cases. The additivity principle, which on 
the other hand can be readily applied to obtain quanti- 
tative predictions, is equivalent within HFT to the hy- 
pothesis that the optimal profiles T(x,t) and j{x,r) 
solution of the variational problem (|TB|) - (|17p are time- 
independent, in which case we recover eq. ((4]) for Q(q). 
In some special cases this approximation breaks down 
for for extreme current fluctuations [3, [TH, [l!| , but even 
so the additivity hypothesis correctly predicts the cur- 
rent LDF in a very large current interval, making it very 
appealing. 



III. THE KMP MODEL 

The system is defined on a ID open lattice with N 
sites [ItJ • Each site models an harmonic oscillator which 
is mechanically uncoupled from its nearest neighbors but 
interact with them through a random process which re- 
distributes energy locally. In this way, a configuration 
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FIG. 1: (Color online) Q(q) for the KMP model as derived 
from the additivity principle, for Tl — 2 and Tr = 1. Notice 
the linear decay for large enough \q\. Vertical lines signal 
the crossover from monotone (\q\ < 7r/3) to non-monotone 
(\q\ > 7r/3) optimal profiles. The Gaussian approximation for 
q w (q), g(q) « -3(|g| - |) 2 / 14 , is a l so shown. 



|q|=3.1 




FIG. 2: (Color online) Optimal T q (x) for different values of \q\, 
both in the monotone and non-monotone regimes, for Tl — 2 
and Tr — 1. The optimal profiles are independent of the sign 
of the current, T q (x) = T- q (x). 



is given by C = {e;, i = 1 . . . N}, where 6 M+ is the 
energy of site z, and the stochastic dynamics proceeds 
through random energy exchanges between randomly- 
chosen nearest neighbors, i.e. (ei,ei+\) — > {e! i ,e' i+1 ) for 
ie [1,N -1] such that 

e i = p(e» + ei+i) 
4+1 = (! -P)(e» + e»+i) > (18) 

with p € [0, 1] a homogeneous random number so + 
e,i + i = e[ + e' i+1 . In addition, boundary sites (i = 1,N) 
may also exchange energy with boundary heat baths at 
temperatures Tl for i — 1 and Tr for i = N , i.e. e^jv — > 
e[ N such that 

e i,N = p(cl,r + ei,iv) (19) 

with e^ij randomly drawn at each step from a 
Gibbs distribution at the corresponding temperature, 
Pk exp(— fik<ik), k = L, R, and p € [0,1] random. For 
T L / T fl KMP proved [l3| that the system reaches a 
nonequilibrium steady state which, in the N — > oo hy- 
drodynamic scaling limit, is described by Fourier's law 
with a nonzero average current 

(«>=-«(T)^M , are [0,1], (20) 
dx 

with k(T) = |, and a linear energy profile 

r Bt (a:)=T i +a;(T fl -T£). (21) 

In addition, convergence to the local Gibbs measure was 
proven in this limit [l7| . meaning that e^, i 6 [l,iV], 
has an exponential distribution with local temperature 
T st [x = i/(N+l)] in the thermodynamic limit. However, 



corrections to Local Equilibrium (LE), though vanishing 
in the iV — > oo limit, become apparent at the fluctua- 
tion level [1^, [53], as we will show below. Moreover, the 
fluctuations of the current in equilibrium (Tl — Tr) are 
described by a(T) = T 2 . It is also worth noticing that 
KMP dynamics obeys the local detailed balance condi- 
tion and is therefore time-reversible [l3j], see Appendix ID] 
In this way we expect the Gallavotti-Cohen symmetry to 
hold in this system, see eq. ([6]). 

The KMP model plays a fundamental role in nonequi- 
librium statistical physics as a benchmark to test new 
theoretical advances, and represents at a coarse-grained 
level a large class of quasi- ID diffusive systems of techno- 
logical and theoretical interest. In this way, understand- 
ing how the energy current fluctuates in the KMP model 
is of central importance to understand current statistics 
in more realistic systems. Furthermore, the KMP model 
is an optimal candidate to test the additivity principle 
because: (i) One can solve eqs. ((4]) and (|5|) to obtain 
explicit predictions for its current LDF, and (ii) its sim- 
ple dynamical rules allow a detailed numerical study of 
current fluctuations. 

In Appendix [X] we apply the additivity formalisms of 
the previous section to study current fluctuations in the 
KMP model. In particular, we use eqs. ((4]) and ([5]) to de- 
rive analytical expressions for the current LDF Q(q) and 
the associated optimal profiles T q (x), see Figs. [TJU In 
this case it can be shown that optimal profiles can be ei- 
ther monotone or non-monotone with a single maximum, 
see Appendix [Al for the explicit calculations. In what fol- 
lows, we compare this set of analytical predictions with 
computer simulation results. 
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FIG. 3: (Color online) Legendre transform of the current 
LDF for the KMP model in one dimension in a temperature 
gradient (top, Tl = 2, Tr = 1) and in equilibrium (bottom, 
Tl = 1.5 = Tr). Symbols correspond to numerical simula- 
tions, full lines to BD theory, and dashed lines to Gaussian 
approximations (see text). Errorbars (with 5 standard devia- 
tions) are always smaller than symbol sizes. The vertical dot- 
ted lines in top panel signal the transition between deviations 
for which the associated temperature profile is monotone (in- 
ner region) or non- monotone (outer region). In equilibrium 
profiles are non-monotone for all current fluctuations. The 
inset in the bottom panel tests the Gallavotti-Cohen relation 
in equilibrium by plotting the difference /i(A) — fj,(— A). 



IV. NUMERICAL TEST OF THE ADDITIVITY 
PRINCIPLE 



The simplicity and versatility of the KMP model al- 
lows us to obtain explicit analytical expressions for Q{q) 
and T q (x) based on the additivity conjecture, see Ap- 
pendix [X] Figs. [T] and [2] show the theoretical current 
LDF and the associated optimal profiles, respectively. 
We find that Pjv(<7, T L , Tr, t) is Gaussian around (q) with 
variance c(T), while non-Gaussian, exponential tails de- 
velop far from (q), with decay rates given by the in- 
verse bath temperatures. Exploring by standard simu- 
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posed. The Gallavotti-Cohen symmetry is satisfied for a wide 
range of A. The inset shows the difference /x(A) — fi( — A — £). 



lations these tails to check BD theory is very difficult, 
since LDFs involve by definition exponentially-unlikely 
rare events. This is corroborated in Appendix [El where 
Q(q) is measured directly but we are unable to gather 
enough statistics in the tails of the current distribution 
to validate or falsify the additivity hypothesis. Recently 
Giardina, Kurchan and Peliti [20] have introduced an effi- 
cient method to measure LDFs in many particle systems, 
based on a modification of the dynamics so that the rare 
events responsible of the large deviation are no longer 
rare [2lj ]. This method yields the Legendre transform of 
the current LDF, fi(X), see eq. (ITUl) . If Ucc is the tran- 
sition rate from configuration C to C of the associated 
stochastic process, the modified dynamics is defined as 
Ucc(X) — Ucc exp(A Jcc), where Jcc is the elemen- 
tary current involved in the transition C — > C . It can be 
then shown (see Appendix |C|) that the natural logarithm 
of the largest ei gen value of matrix U(X) gives ^t(A). The 
method of Ref. [20J thus provides a way to measure fi(X) 
by evolving many copies or clones of the system using 
the modified dynamics U (A) , see Appendix [C] 

We applied the method of Giardina et al. to mea- 
sure (j,(X) for the ID KMP model with TV = 50, T L = 2 
and Tr = 1, see Fig. [3J top panel. The agreement 
with BD theory is excellent for a wide A-interval, say 
—0.8 < A < 0.45, which corresponds to a very large 
range of current fluctuations, see inset to Fig. [TT] in 
Appendix [Cj Moreover, the deviations observed for ex- 
treme current fluctuations are due to known limitations 
of the algorithm [l4], l20T - |22l ]. so no violations of addi- 
tivity are observed. Notice that the spurious differences 
seem to occur earlier for currents against the gradient, 
i.e. A < 0. In fact, we can use the Gallavotti-Cohen 
symmetry, which in A-space now reads /x(A) = A — £) 
with £ = (T^ 1 - T^ 1 ), to bound the range of validity 
of the algorithm: Violations of the fluctuation relation 
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indicate a systematic bias in the estimations provided by 
the method of Ref. [H, see also (H). Fig. U shows that 
the Gallavotti-Cohen symmetry holds in the large cur- 
rent interval for which the additivity principle predictions 
agree with measurements, thus confirming its validity in 
this range. However, we cannot discard the possibility 
of an additivity breakdown for extreme current fluctua- 
tions due to the onset of time-dependent optimal profiles 
expected in general in HFT 4], although we stress that 
such scenario is not observed here. 

We also measured the current LDF in canonical equi- 
librium, i.e. for Tl = Tr = f.5, see the bottom panel 
in Fig. [3] The agreement with BD theory is again ex- 
cellent within the range of validity of our measurements, 
which expands a wide current interval, see inset to Fig. [3J 
and the fluctuation relation is verified except for extreme 
currents deviations, where the algorithm fails to provide 
reliable results. Notice that, both in the presence of a 
temperature gradient and in canonical equilibrium, /i(A) 
is parabolic around A = meaning that current fluctu- 
ations are Gaussian for q « (q), as demanded by the 
central limit theorem, see eqs. (|A15|) - (|A16|) in Appendix 
lAl This observation is particularly interesting in equi- 
librium, where canonical and microcanonical ensembles 
behave differently (see below). 

The additivity principle leads to the minimization of 
a functional of the temperature profile, T q (x), see eqs. 
Q and ([5]). A relevant question is whether this opti- 
mal profile is actually observable. We naturally define 
T q (x) as the average energy profile adopted by the sys- 
tem during a large deviation event of (long) duration t 
and time-integrated current qt, measured at an interme- 
diate time 1 <C t <C t, i.e. T q (x) = T q nid (x) . Fig. 
[5] shows the measured T™ ld (x) for both the equilibrium 
and nonequilibrium settings, and the agreement with BD 
predictions is again very good in all cases, with discrep- 
ancies appearing only for extreme current fluctuations, 
as otherwise expected. See also Fig. [TJ]in Appendix iBl 
This confirms the idea that the system indeed modifies 
its temperature profile to facilitate the deviation of the 
current, validating the additivity principle as a power- 
ful conjecture to compute both the current LDF and the 
associated optimal profiles. Our numerical results show 
also that optimal profiles are indeed independent of the 
sign of the current, T\(x) = T-\-s(x) or equivalently 
T q (x) = T-q(x), a counter- intuitive symmetry resulting 
from the reversibility of microscopic dynamics. Notice 
that in the equilibrium case (Tl = Tr) optimal tem- 
perature profiles are always non-monotone with a single 
maximum for any current fluctuation q ^ (q) (the sta- 
tionary profile is obviously flat). This is in stark con- 
trast to the behavior predicted for current fluctuations 
in microcanonical equilibrium, i.e. for a one-dimensional 
closed diffusive system on a ring 0, EE EH • In this case 
the optimal profiles remain fiat and current fluctuations 
are Gaussian up to a critical current value, at which 
profiles become time-dependent (traveling waves) [l9| . 
Hence current statistics can differ considerably depend- 
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FIG. 5: (Color online) Excess temperature profiles for differ- 
ent current fluctuations (O), f° r a system subject to a tem- 
perature gradient (top, Tl = 2, Tr = 1) and in equilibrium 
(bottom, T L = 1.5 = Tr). In all CELS6S, 3. greement with BD 
theoretical predictions (lines) is very good within the range of 
validity of the computational method. Dotted symbols cor- 
respond to midtime profiles obtained from endtime statistics 
(see text). 



ing on the particular equilibrium ensemble at hand, de- 
spite their equivalence for average quantities in the ther- 
modynamic limit. Finally, notice also that equilibrium 
optimal profiles are symmetric with respect to x = 1/2, 
as expected since Tl = Tr. 

For small enough current fluctuations around the av- 
erage, q (q) with (q) = 1/2 for T L = 2 and T R = 1, 
BD theory predicts the limiting behavior 



T q (x) - T st (x) 
2q-l 



1 , 
7^(1 



x)(5 



0(2g-l). (22) 



Fig. |5] confirms this scaling for T q (x) and many different 
small current fluctuations around the average. In partic- 
ular, it shows data obtained both from standard simula- 
tions (see Appendix [B| and using the advanced method 
of Ref. [H. 



x=i/(N+1) 



FIG. 6: (Color online) Scaling plot of the excess profiles for 
small current fluctuations. Here we plot results obtained from 
standard simulations (solid circles) and the advanced algo- 
rithm of Ref. [20| (open squares), as well as the theoretical 
prediction (line). 

It is also interesting to study the statistics of config- 
urations both during a large deviation event and at the 
end. They differ due to final transient effects which de- 
cay exponentially fast, but a connection exists between 
both regimes which highlights the symmetry of midtime 
statistics resulting from the reversibility of microscopic 
dynamics (a symmetry akin to the fluctuation relation). 
Reversibility in stochastic dynamics stems from the con- 
dition of local detailed balance [13|, which implies a re- 
lation between the forward modified dynamics for a cur- 
rent fluctuation, U(X), and the time-reversed modified 
dynamics for the negative fluctuation, U T (—X — £), see 
eq. (|D6|) in Appendix [D] This can be used to derive 
a relation between midtime and endtime statistics (see 
Appendix |D|) . 

pend^pend ! C \ 
pnnd {c) = A A K ) -X- £ \ ) 23 

Pc 

Here P% nd (C) [resp. P^ id (C)} is the probability of con- 
figuration C at the end (resp. at intermediate times) 
of a large deviation event with current-conjugate pa- 
rameter A, and p°Q = exp[— J3i=i A e *] is an effective 
weight for configuration C — {ei,i = 1...N}, with 
Pi = T^ 1 + £ ir\ , while A is a normalization con- 
stant. Eq. (|23|) implies that configurations with a sig- 
nificant contribution to the average profile at interme- 
diate times are those with an important probabilistic 
weight at the end of both the large deviation event and 
its time-reversed process. An important consequence of 
eq. K23J) is hence that P£ id (C) = P mi x d _ £ {C), or equiva- 
lently P g mid (C) = P™ d (C), so midtime statistics does 
not depend on the sign of the current. This implies 
in particular that T q nid (x) = T™ d (x), but also that 
all higher-order profiles (e n (x)) q and spatial correlations 



FIG. 7: (Color online) Semilog plot of local energy histograms 
along the chain for different values of A, at the end of the large 
deviation event. Notice that, in all cases, energy distributions 
are very close to exponential. 

(e n (xi) . . . e n (x m )) q are independent of the current sign 
Vn, m. 

The above connection allows us to relate midtime and 
endtime profiles for a given current fluctuation. For that 
we need additionally a local equilibrium (LE) hypothe- 
sis, i.e. we now assume that spatial correlations at the 
end of a large deviation event are weak enough so the 
distribution P™ d (C) can be approximately factorized, 
P° nd (C) « n^ 1 P™ d (e i ). In this way we obtain a lo- 
cal equilibrium picture with local temperature parame- 
ter T^ nd (x = jyVj). This hypothesis can be numerically 
justified by measuring, at the end of the large deviation 
event, local energy distributions along the chain for dif- 
ferent values of A, see Fig. [7] In all cases the distribution 
is compatible with local equilibrium to a large degree of 
accuracy. Using eq. (|2"B")1 and the LE hypothesis we thus 
find 

Tmid(x) = Tr d (x)Tlf_ £ (x) 

x { T A ™ d (x) + T™ d _ f (x) - p x T™ d (x) T- d _ £ (x) ■ 

(24) 

Fig. [8] shows endtime profiles T^ nd (x) measured both 
in equilibrium (bottom) and nonequilibrium (top) con- 
ditions for different values of A. These profiles are 
clearly asymmetric upon current inversion, T™ d (x) ^ 
TZ nd _ £ (x), and most interestingly they show boundary 
resistance which depends on A and on the particular def- 
inition for the elementary current, see 22] . In the equilib- 
rium case the symmetry T° nd (x) — T™ d (f — x) resulting 
from the reflection invariance in this case (Tl = Tr) is 
apparent in Fig. [8] (bottom). Fig. [5] also shows midtime 
profiles obtained from the measured T^ nd (x) via eq. (|24|) . 
The agreement with theoretical predictions and direct 
measurements of midtime profiles is good, though dis- 
crepancies appear for large enough current fluctuations, 
pointing out that corrections to LE are weak but increase 
for large current deviations. We show below that these 
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FIG. 8: (Color online) Excess temperature profiles measured 
at the end of the large deviation event for different values 
of A, both in the presence of a temperature gradient (top), 
Tl = 2 and Tr = 1, and in canonical equilibrium (bottom), 
Tl = 1-5 = T R . Notice that in all cases T^ nd {x) ^ Tlf_ £ (x), 
although for the equilibrium case the symmetry Tl nd (x) = 



Tl n x d (l 



x) is apparent. 



corrections are also present for small current fluctuations 
and can be measured. 

We can now explore physics beyond the additivity con- 
jecture by studying fluctuations of the system total en- 
ergy, e(C) = ./V -1 for which current theoretical 
approaches cannot offer any general prediction. An ex- 
act result by Bertini, Gabrielli and Lebowitz (BGL) [27} 
predicts that 



m 2 (e) 



LE 



(25) 



where m^e) — N((e 2 ) — (e) 2 ) is the variance of the total 
energy in the nonequilibrium steady state (NESS), m^ E 
is the variance assuming a local equilibrium (LE) prod- 
uct measure, and the last term reflects the correction to 




FIG. 9: (Color online) Fluctuations of the total energy per 
site versus A for both equilibrium (□, Tl = 1.5 = Tr) and 
nonequilibrium (O, Tl, = 2, Tr = 1) conditions. The lines 
stand for predictions based on the additivity principle plus 
a local equilibrium hypothesis. Inset: Average energy per 
site and BD prediction in both situations. Notice that, as 
before, deviations observed in all cases for extreme current 
fluctuations are spurious and result from known limitations 
of the method of Ref. EH. 



LE due to weak long-range correlations in the NESS [27J , 
which in this case results in the enhancement of energy 
fluctuations. Corrections to LE vanish in the thermo- 
dynamic limit but extend over macroscopic distances (of 
order N) , giving rise in general to a non-local current 
LDF [27]. In our case, 



LE 



1 



T L T R 



Tl) = 



7 



2.3333, (26) 



29/12 » 2.4166. Fig. [9] plots m 2 (e,A) = 
21 as a function of A for both equilibrium and 



while 77i2 
N[{e*) x - 

nonequilibrium conditions, showing a non-trivial, inter 
esting structure which both BD theory and HFT cannot 
explain. One might obtain a theoretical prediction for 
TO2(e, A) by supplementing the additivity principle with 
a LE hypothesis, 



Px{C) ocn^exp 

which results in 

ml E (e,X)-- 



T\(x) 2 dx . 



(27) 



(28) 



This prediction agrees qualitatively with the observed be- 
havior, though fine quantitative differences are apparent, 
see Fig. [9l as otherwise expected. In particular we find 
that, out of equilibrium, m^ E {e 1 0) 2.33 as corresponds 
to a LE picture, and in contrast to the measured value 
m<x{e, 0) = 2.422(14) in Fig. [SJ which compares nicely 
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with the exact BGL result 29/12 (recall that A = cor- 
responds to q = (q)). This shows that, even though LE 
is a sound numerical hypothesis to obtain T\(x) from 
cndtime statistics for small and moderate current fluc- 
tuations, see Fig. [5] and eq. (|24l) . corrections to LE 
become apparent at the fluctuating level even for small 
current fluctuations. This is also shown in Fig. [15] in 
Appendix [B] where fluctuations of the total energy un- 
der nonequilibrium conditions are studied in standard 
simulations. On the other hand, in the canonical equilib- 
rium case (T^ = 1.5 = Tr) no corrections to LE show up 
for A = (i.e., for q = (q) =0), as expected. However, 
as soon as q ^ (q), deviations of m 2 (e,A) from the LE 
prediction m^fe, A) are observed, thus showing that lo- 
cal equilibrium is broken at the fluctuating level even for 
equal bath temperatures. 

Finally, the inset to Fig. |H] shows the average energy 
per site as a function of A, together with the prediction 
based on the additivity principle, (e)\ — J Q T\(x)dx. 
Agreement is again very good in the large range of cur- 
rents explored. It is interesting to note that in or- 
der to sustain a current fluctuation above the average, 
q > (q) or equivalently A > 0, the nonequilibrium sys- 
tem (Tj, > Tr) has always a larger average energy than 
its equilibrium counterpart (Tl = Tr), while the reverse 
holds for current fluctuations below the average, q < (q), 
see inset to Fig. [S] 



V. JOINT FLUCTUATIONS OF THE CURRENT 
AND THE PROFILE 

For long but finite times, the profile associated to a 
given current fluctuation is subject to fluctuations itself. 
These joint fluctuations of the current and the profile are 
again not described by the additivity principle, but we 
may study them by extending the additivity conjecture. 
In this way, we now assume that the probability to find 
a time-integrated current q/N and a temperature profile 
T q (x) after averaging for a long but finite time t can be 
written as 
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FIG. 10: (Color online) Finite-time profile fluctuations 
{Wq(x)). Blue squares are the numerical evaluation of the se- 
ries expansion from the extended BD theory (see text). Black 
and red circles are standard simulation results for g's in the 
interval shown in the figure, N = 50 and t — 4000 and t — 10 6 . 



For large enough t, the joint probability of q and T) q (x) 
can be written as 



W N [q,r) q (x);t] 



exp 



1 



dxdyA q (x,y)r] q (x)T]g(y) 



(32) 

where PAr(g,i) is defined in equation ([TJ, together with 
eqs. (j4j and ([5]). The integral kernel is 



1 dT q d 



1 d 2 



2T 9 3 dx dx AT 2 dx 2 



- 2 



K(q)q 2 



rp 2 



S(x - y) . 



(33) 



One can show that the kernel A q {x, y) is symmetric with 
respect to x and y. In order to check the above joint 
probability distribution, we studied the observable 



(f 2 (x)) - T 2 (x) = ( V 2 (x)) = A-\x,x) , (34) 



where 



W N [j-,T 9 (x);t] 



exp 



t 

~N 



G[q,T q (x)\ 



(29) 



A q 1 (%> y) = ^2<t> n 1 Vn(x; q)v n (y; q) , (35) 



where now 



G[q,T q (x)] 



1 [q + Klf q (x)]f>(x)Y 
2a[f q (x)} 



dx . 



(30) 



Notice that here no minimization with respect to temper- 
ature profiles is performed, see eq. @. In this scheme 
the profile obeying eq. ([5]), i.e. the one which minimizes 
the functional Q, is the classical profile Tq(x). For a given 
q value we can make a perturbation of T q (x) around its 
classical value, 



T q (x) = T q {x) + r) g (x) . 



(31) 



and v n {x] q) and 4>n are the eigenvectors and eigenvalues 
of kernel A q , respectively, 



dxA q (x, y)v n (x;q) = <p n v n {y; q) , 



(36) 



with v n (0; q) = = v n (l; q). For q = (q) = 1/2 (nonequi- 
librium conditions, Tj, = 2, Tr = 1) we were able to solve 
the eigenvalue equation, yielding 



v n (x;l/2) = BT 1/2 (x) 3 / 2 ^J_ 3/ ^ n T 2 )J 3/ SnT 1/2 (x) 2 } 
- J3/4(0„r|)J_ 3 /4[^T 1/2 (x) 2 ]l, (37) 
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where <f) n = (<j) n N /t) 1 / 2 / (T L - T R ), J's are the Bessel 
functions and B is the normalization factor that is ob- 
tained by requiring 



dxv n (x;l/2) 2 = 1 



Finally, (f> n are the solutions of the equation 



(38) 



<- / 3/4(</ > n7 1 !)J_3 /4 ((/> n T^) = J_ 3/4 (4> n Tl)J 3/4 ((j) n T^) . 

We compare in Fig. [TU] the numerical evaluation of 
A~h(x,x) (where we have computed 10, 15, 30, 50, 100 
and 200 terms of the series and extrapolated to n —> oo) 
with the standard simulation results for N — 50 and 
t = 4000 and t = 10 6 . We observe a good agreement 
between theoretical and simulation results. Notice that 
we average over a small q- window around q — 1/2 in 
simulations. These results show that the BD functional 
G[q,T q (x)] of eq. (I3TJ1) contains the essential information 
on the joint fluctuations of the current and the average 
profile, extending the validity of the additivity principle 
to finite-time situations. 



VI. CONCLUSIONS 

In this paper we have confirmed via extensive com- 
puter simulations the validity of the additivity princi- 
ple for current fluctuations in the ID Kipnis-Marchioro- 
Pressuti model of energy transport. In particular, we 
found that the current distribution shows a Gaussian 
regime for small current fluctuations and non-Gaussian, 
exponential tails for large deviations of the current, such 
that in all cases the fluctuation relation holds. We ver- 
ified the existence of a well-defined temperature profile 
associated to a given current fluctuation, different from 
the steady-state profile and invariant under current rever- 
sal. In addition, we extended the additivity conjecture 
to joint current-profile fluctuations. 

Our results thus strongly support the additivity hy- 
pothesis as an important tool to understand current 
statistics in diffusive systems, opening the door to a gen- 
eral approach to a large class of nonequilibrium phenom- 
ena based on few simple principles. Our confirmation 
does not discard however the possible breakdown of ad- 
ditivity for extreme current fluctuations due to the onset 
of time-dependent profiles, although we stress that this 
scenario is not observed here and would affect only the far 
tails of the current distribution. In this respect it would 
be interesting to study the KMP model on a ring, for 
which a dynamicphase transition to time-dependent pro- 
files is expected @,[ll,[l{|. Also interesting is the possible 
extension of the additivity principle to low-dimensional 
systems with anomalous, non-diffusive transport proper- 
ties [lj| , or to systems with several conserved fields or in 
higher dimensions. 



Appendix A: Predictions using the Additivity 
Principle 

In this appendix we use the KMP model values for 
k(T) = \ and cr(T) = T 2 in eqs. Q and © to de- 
rive explicit predictions for the current large deviation 
function in this model and the associated optimal tem- 
perature profiles. In what follows we assume Tl > Tr 
without loss of generality. The differential equation for 
the optimal profile in the KMP case reads 



dT q (x) 
dx 



4q 2 {l + 2T 2 (x)K(q 2 )} 



(Al) 



Here two different scenarios appear. On one hand, for 
large enough K(q 2 ) the rhs of eq. (|A1|) does not vanish 
V.t G [0, 1] and the resulting profile is monotone. In this 
case, the optimal profile obeys 



dT q (x) 
dx 



2\q\Jl + 2T 2 (x)K(q 2 ) 



(A2) 



On the other hand, for K{q 2 ) < the rhs of eq. (IA1|) may 
vanish at some points, resulting in a T q (x) that is non- 
monotone and takes an unique value T* = \f— \/2K(q 2 ) 
in the extrema. Notice that the rhs of the above equation 
may be written in this case as 4g 2 [l — (T q (x)/T*) 2 ]. It 
is then clear that, if non- monotone, the profile T q (x) can 
only have a single maximum T q (x*) — T* because: (i) 
Tq{x) < T* Vx £ [0, 1] for the profile to be a real function, 
and (ii) several maxima are not possible because they 
should be separated by a minimum, which is not allowed 
because of (i). In this case 



dT q (x) 
dx 



-2M 



T q (x) 



-2\q\ K l 



T Q (x) 



x < x 



X > X 



(A3) 



This leaves us with two separated regimes for current 
fluctuations, with the crossover happening for \q\ = 

This crossover current may be ob- 



2 



"'(ft) 



tained from eq. (|A13|) below by letting T* -> T L 



1. Region I: \q\ < ^ [f - sin^ 1 



In this region the optimal profile T q (x) is monotone in 
x S [0, 1]. Eq. (HJ then leads to 



q ( 1 1 



2 V T R Tt 



q 2 K(q 2 ' 



(A4) 



|g| / y/1 + 2K(qi)T 2 y/1 + 2K(g*)T 2 

2 Tr Tfj 
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where K(q 2 ) is a constant defined by the boundary con- 
ditions. The optimal temperature profile T q (x) in this 
regime is the solution of the following implicit equation 



2.r|g| = 



1 



In 



T q (x) + ^T q (x) 2 + 2K(q'i) 



whenever K(q 2 ) > 0, or rather 



(A5) 



sin 



2x\q\ = 



y/-2K(q 2 )T L 



y/-2K(q 2 )T q (x) 



y/-2K(q*) 



(A6) 

in the case — 157^ < K(q 2 ) < 0, see eq. (|A2j) . Making 

x = 1 and T q (x = 1) = Tr here, we obtain the implicit 
equation for the constant K(q 2 ). 

Some times it is interesting to work with the Legen- 
dre transform of the large deviation function, /i(A) — 
N^ 1 max, [G(q) + Xq] = Q(q ) + Xq a , with q a (X) given by 
d q G(q ) + A = 0, and where now -T^ 1 < A < T£ l . It 
then follows 



MA) = -^bo(A)] 2 



where 



2|?o(A)| = 



In 



V2K{X) 
when K(X) > 0, or instead 



T 



R + \J T R + JElTj 



(A7) 



(A8) 



2|g (A)| 



y/-2K(X)T L 



y/-2K(X)T R 



y/-2K{\) 



(A9) 



in the case —t^x < K(X) < 0, and the constant K(X) = 
K[q (X) 2 } is solution of the implicit equation 



i i-i 
2 \T R T L 



(A10) 



sgn [q Q (X)} 



y/1 + 2K(X)T 2 R y/1 + 2K(X)T 1 



Tr 



T l 



The optimal profile for a given A is just T\(x) = 
Tq o (v\(x). In A-space, monotone profiles are expected 
for A G [A_ , A+] where A± = -{Tr 1 - T^)/2 ± 

y/1 - (Tr/T L ) 2 /(2Tr). 



2. Region II: \q\ > [f - sin" 1 (^)] 

In this case the optimal profile is non-monotone with 
a single maximum T* = T q (x*), see eq. (|A3[) . In this 



regime K(q 2 ) < 0, and T* = l/y/-2K(q 2 ). It follows 



7r — sin 



— sm 



(All) 



2 Vr* 



1 



Tr, 



Tr 



The optimal profile solution of eq. (|A3p is given by 



q 

2kl 



-1 



1 + 



q 



sm 



— sm 



± q 



X < X 



(A12) 



a; > x 



At the location of the profile maximum, x = x* , both 
branches in the above equation must coincide and this 
condition provides equations for both x* and T* 



sm 

2 



Tl 



7T — sin 



Tl 



(A13) 



(A14) 



As in Regime I, we find for the Legendre transform 
/i(A) = -N- 1 K{X)q {X) 2 = (2N)- 1 [q (X)/T*} 2 , with 
g (A) defined in eq. (|A"l3|) , T* = T* (A) , and A given 
as in eq. (IA10[) but with the notation change K(X) — > 
— 1/[2(T£) 2 ]. Non- monotone profiles are then expected 
for A G {—T^ 1 , A_) U (A+, T^ 1 ]. 

Fig. [T]in the main text shows the predicted Q (q) for 
the KMP model. Notice that the large deviation function 
is zero for q = (q) = (Tl —Tr) /2, and negative elsewhere. 
Moreover, for large current fluctuations it decays linearly, 
Q(q) — > —q/TR.L for \q\ 3> (q). For a small positive 
current fluctuation, K(q 2 ) — > and 



G(q) « - 
which translates into 



Tl - T R 



2(T 2 + T l Tr + T 2 /) 



(A15) 



X 

2N 



(T L - T R ) 



(Tl 



TlTr + T 2 R ) 



(A16) 

for the Legendre transform. Therefore the probability 
of small current fluctuations is Gaussian in q while it 
becomes exponential for large enough deviations from 
the average, see eq. ([1]). It is easy to show that the 
Gallavotti-Cohen symmetry holds, with 



Q(q) - G(-q) = 2q - 



K(T) 

(T) 



1 

T^. 



Tr 



, (A17) 
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FIG. 11: (Color online) Top panel: Constant K as a function 
of q for Tr = 2 and Tl = 1. Bottom panel: the same constant 
as a function of A. The inset shows the current q conjugated 
to A. 



or equivalently 



m (A) = m(-A-£), 



(A18) 



with E = (Tft 1 — T7 l ). Fig. din the main text shows the 
optimal temperature profiles for different current devia- 
tions. Notice that the optimal profile is independent of 
the sign of the current, i.e. T q {x) = T- q (x), reflec ting th e 
time-reversal symmetry of microscopic dynamics [121 1 1 31 ) . 
Finally, Fig. 1111 shows, for information purposes, the in- 
tegration constant K as a function of both q and A, as 
well as A-dependence of g (A). 



Appendix B: Standard Simulations 

In order to see how far standard simulations can go in 
evaluating current large fluctuations, and to cross-check 
our results with the more advanced simulation methods 
described in Appendix [C] we performed a large num- 
ber of steady-state simulations of long duration £, with 
Tl, = 2 and Tr = 1, measuring the total time-integrated 
current Q t = qt and accumulating statistics for q. Fig. 
[T2"l shows the measured G{q) obtained for different system 




FIG. 12: (Color online) Q(q) measured for different system 
sizes iV and measurement times t (see text), with Tl = 2 
and Tr — 1 fixed. Lines correspond to BD theory and the 
Gaussian approximation. 



sizes iV and durations t. Our simulations for N = 1000 
and different times t < N 2 follow closely the Gaussian 
law Q(q) ~ — 3(q — l/2) 2 /14 obtained from the first two 
moments prescribed by the additivity principle in this 
case, namely 



mi 



m 2 



T L - Tu 



Tl + T L T R + Tj 



This Gaussian behavior is expected for small fluctuations 
around the average current, see eq. (|A15[) . but deviations 
away from Gaussianity should be already observed in the 
current range studied, see the theoretical prediction. In 
particular, the theoretical Q{q) implies a nonzero third 
central moment, but we have not found numerical evi- 
dence of such a deviation for N — 1000. This lack of 
structure stems from the relatively short duration of the 
simulations for N = 1000, i.e. our results are not in the 
diffusive regime (t < N 2 here) and therefore we have not 
reached the asymptotic behavior. 

We performed two set of simulations in the diffusive 
regime t > N 2 , namely TV = 50 with t = 10 6 and 
t = 4000. In the first case there were no events out- 
side the current interval q € [0.45,0.56], for which the 
BD prediction is numerically indistinguishable from the 
Gaussian one. On the other hand, the case N = 50 and 
t = 4000 shows systematic deviations from Gaussian be- 
havior, seemingly compatible with BD theory, see Fig. 
1121 However, large errorbars resulting from the difficulty 
of gathering statistics in this rare-fluctuation regime do 
not allow us to exclude Gaussian behavior. In this way, 
standard simulation results are inconclusive, as otherwise 
expected, and the more refined simulation techniques of 
Appendix [Cl are called for. 

We also tested the Gallavotti-Cohen relation in stan- 
dard simulations for our system. This symmetry implies 



13 



0,0025 



+t. 0,002 
cr 



Q_ 0,0015 

-I— ' 

cr 



0,001 



0^ 

£= 



r- — 0,0005 



1/2 (Fluctuation Theorem) 




Linear regression: 

1-tOOO, y=0.456(0.001): 
t=2000, y=0.470(0.001): 
t=3000, y=0.472(0.005): 
t=4000, y=0.470(0.006): 



FIG. 13: (Color online) Test of the fluctuation theorem of 
Gallavotti and Cohen. Here we explore N = 50 and different 
maximum times t. If BD theory holds a slope 1/2 is expected, 
while Gaussian behavior involves a slope 3/7. 
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FIG. 14: (Color online) Excess average profiles during a large 
deviation event for small current fluctuations, as measured in 
standard simulations. Agreement with BD theoretical predic- 
tions (lines) is excellent. 



that 



t^oct P N (-q,T L ,T R ,t) 



= £q, 



(Bl) 



where £ = (T^ 1 - T^ 1 ) = 1/2 in this case. Notice that 
if we assume Ppf(q,TL,Tfi;t) to be Gaussian with the 
moments defined above, then one expects £ = 3/7. Fig. 
[T31 shows the above quotient as measured for N = 50 
and different values of t. It shows a systematic deviation 
from Gaussian behavior which increases with t. However, 
we do not see clearly £ — 1/2, and this means again 
that our standard simulations are still far from the true 
asymptotic regime in t. 

Another prediction of the additivity principle concerns 
the existence of an optimal temperature profile that the 
system adopts in order to facilitate a given current fluctu- 
ation. We measured in standard simulations the average 
energy profile during a current large deviation event, ob- 
taining the results shown in Fig. [14] As above, only for 
small current fluctuations we could gather enough statis- 
tics for the data to be significative. In any case, the the- 
oretical optimal profiles compare nicely with data, con- 
firming the existence of a well-defined temperature profile 
for each current deviation. 

We also measured the fluctuations of the total energy 
in standard simulations. Fig. [T^lshows our results in this 
case. In particular, we measured 777.2(e) = 2.4 (1) for N = 
50 and a maximum time t = 4000 and 777.2(e) = 2.42 (2) 
for t = 10 6 , in agreement with eq. (|26|) . This figure also 
shows 777,2 (e, l) and E (e, q) build from simulation data 
for T q {x). As in Fig. [5J we see a clear deviation from local 
equilibrium and a well defined structure not predicted 
by BD theory. Notice that, again, values of 7772 (e,q) for 
q = 1/2 coincide with the expected average values with 
no current constraint. The data shown in this figure agree 
nicely with those measured with the advanced technique 
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FIG. 15: (Color online) Fluctuations of total energy vs q mea- 
sured in standard simulations for t = 10 6 (O) an d LE results 
(□). Inset: Similar results for t = 4000. Notice the non-trivial 
structure. 



in the studied range, see Fig. [SI 



Appendix C: Evaluation of Large-Deviation 
Functions 



Large deviation functions are very hard to measure in 
experiments or simulations because they involve by defi- 
nition exponentially- unlikely events, see eq.([T]). Recently, 
Giardina, Kurchan and Peliti [20| have introduced an ef- 
ficient algorithm to measure the probability of a large 
deviation for observables such as the current or density 
in stochastic many-particle systems. The algorithm is 
based on a modification of the underlying stochastic dy- 
namics so that the rare events responsible of the large de- 
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viation are no longer rare, and it has been extended for 
systems with continuous-time stochastic dynamics [2lj |. 
Let Ucc be the transition rate from configuration C to 
C . The probability of measuring a time- integrated cur- 
rent Q t after a time t starting from a configuration Cq 
can be written as 

t-i 

P(Qt,t;C a )= £ U CtOt - 1 ..Uc 1 c S(Q t -Y l Jc„ 1 c l .), 

CfCi k=0 

. ( C1 ) 

where Jcc is the elementary current involved in the 
transition C — > C For long times we expect the infor- 
mation on the initial state Co to be lost, P(Qt,t; Cq) — > 
P(Qt,t). In this limit P(Q t ,t) obeys the usual large 
deviation principle P{Q t ,t) ~ exp[+tJ 7 (q = Qt/t)]. In 
most cases it is convenient to work with the moment- 
generating function of the above distribution 

U(X,t) = ]Tc A *P(Q t ,t) (C2) 
Qt 

Ct-.d 

For long t, we have II(A,t) — > exp[+t/i(A)], with jit(A) = 
maxq [J-(q) + Xq] ■ We can now define a modified dynam- 
ics, Ucc = e XJ °' c Ucc, so 

n(A,t)= ]T u Ct c t ^...u Cl c . (C3) 

C t ...Ci 

This dynamics is however not normalized, J^c ^cc 7^ 1- 
We now introduce Dirac's bra and ket notation, use- 
ful in the context of the quantum Hamiltonian formalism 
for the master equation [H [H[, see also [U H3- The 
idea is to assign to each system configuration C a vector 
C) in phase space, which together with its transposed 
vector (C|, form an orthogonal basis of a complex space 
and its dual [23|, HH. For instance, in the simpler case 
of systems with a finite number of available configura- 
tions (which is not the case for the KMP model), one 
could write \C) T = (C\ = (. . . . . . 0, 1, . . . . . .), i-e. 
all components equal to zero except for the component 
corresponding to configuration C, which is 1. In this no- 
tation, Ucc = (C'\U\C), and a probability distribution 
can be written as a probability vector 

\P(t))=Y / p (C,t)\C), 

c 

where P(C,t) — (C\P(t)) with the scalar product 
{C'\C) — Sec- K ( s l = (1 • • • 1): normalization then 
implies (s\P(t)) = 1. 

With the above notation, we can write the spectral 
decomposition U(X) = £V . e A ^( A) \Af(X))(Af(X)\, where 
we assume that a complete biorthogonal basis of right 
and left eigenvectors for matrix U exists, U\A^(X)) = 

e A iW|Af(A)) and (Aj(X)\U = e A ^^(Af(X)\. Denoting 
as e A ( A ) the largest eigenvalue of U(X), with associated 




Mj 



FIG. 16: (Color online) Sketch of the evolution and cloning of 
the copies during the evaluation of the large deviation func- 
tion. 

right and left eigenvectors |A fl (A)} and (A L (A)|, respec- 
tively, and writing Il(X,t) = ^2 Ct (C t \U t \C ), we find for 
long times 

U(X,t)^c+ tA W(A L (X)\C ) (^T(C t \A R (X))^ . 

(C4) 

In this way we have fi(X) = A(A), so the Legendre trans- 
form of the current LDF is given by the natural logarithm 
of the largest eigenvalue of U(X). In order to evaluate this 
eigenvalue, and given that dynamics U is not normalized, 
we introduce the exit rates Yq = Y^c Ucc, an d define 
the normalized dynamics U' c , c = Yq Ucc ■ Now 

U(X,t)= £ Y Ct _ x U' CtCt _ 1 ...Y Co U' CxCo (C5) 
C*...Ci 

This sum over paths can be realized by considering an 
ensemble of M > 1 copies (or clones) of the system, 
evolving sequentially according to the following Monte 
Carlo scheme [20(: 

I Each copy evolves independently according to mod- 
ified normalized dynamics U' c , c . 

II Each copy m € [1,M] (in configuration Ct[m] at 
time t) is cloned with rate Y Ct r OT i . This means 
that, for each copy m 6 [1, M], we generate a num- 
ber Kc t [ m ] — [Yc t [m]\ + 1 of identical clones with 
probability Y Ct[m] - [ic t [m]J ; or K c t [m\ = ITc«[m]J 
otherwise (here \x\ represents the integer part of 
x). Note that if K Ct [ m ] — the copy may be killed 
and leave no offspring. This procedure gives rise to 
a total of M[ — Yl m =i ^c t [m] copies after cloning 
all of the original M copies. 

Ill Once all copies evolve and clone, the total num- 
ber of copies M' t is sent back to M by an uniform 
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cloning probability X t = MjM[. 

Fig. [TH] sketches this procedure. It then can be shown 
that, for long times, we recover /z(A) via 



MA) 



for t > 1 



(C6) 



To derive this expression, first consider the cloning dy- 
namics above, but without keeping the total number of 
clones constant, i.e. forgetting about step III. In this 
case, for a given history {C t , Ct-i . . . C\, Co}, the num- 
ber M(Ct . ■ ■ Cq, t) of copies in configuration Ct at time t 
obeys Af(C t ...C ,t) = Y Ct _ 1 U' Ct c t .J^{O t -i...C Q ,t - 
1), so that 



Af(C t ...C ,t) = Yc^U'c^r ■ ■■Y Ca U' CiCa N{C 0l Q) . 

(C7) 

Summing over all histories of duration t, see eq. (|C5I) . 
we find that the average of the total number of clones 
at long times shows exponential behavior, (N(t)) — 
Ec t ...c 1 N(C t ...C ,t) ~ Af(Co,0)exp[+i/i(A)]. Now, 
going back to step III above, when the fixed number 
of copies M is large enough, we have X t = (Af(t — 
1))/ (Af(t)) for the global cloning factors, so X t ■■ ■ Xi = 
Af(Co, 0) / (Af(t)) and we recover expression (|C6[) for /i(A). 

In this paper we used the above method to measure the 
current LDF for the Kipnis-Marchioro-Presutti model in 
one dimension, described in Section ITTTl For this model 
the transition rate from a configuration C — {e\ . . . e./v} 
to another configuration C' y — {ei . . . e' y ,e' y+1 . . . cn}, 
with y £ [0, N] and the pair (e' , e' y+1 ) defined as in eqs. 
(IT51) - (fTTjl) . can be written as 



y(=[l,N-l] 



/3_e^- ei , 

u C 'c = { n + i El ^- max ( e i' e i)] ' y = ° 
N + i El [^+ max ( e ^' e jv)] > v = N - 

Here Ei(x) = — Ei(— x), where Ei(x) is the exponential 
integral function, or 



Ei (X) 



du- 



(C8) 



It appears when integrating over all possible pairs 
(p,^l,r) that can result on a given e' 1N , respectively, 
see eq. ([TUf in Section Unl It is easy to show that Uc y c 
is normalized as it should, so V„, Ucc = 1- 

y y 

In order to measure current fluctuations we need to 
provide a microscopic definition of the energy current in- 
volved in an elementary move. There are many different 
ways to define this current: the energy exchanged per 
unit time with one of the boundary heat baths, the cur- 
rent flowing between two given nearest neighbors, or its 
spatial average, etc. Assuming that ener gy c annot ac- 
cumulate in the system ad infinitum^, M. |25|. all these 



definitions give equivalent results for the current large 
deviation function in the long time limit. However, this 
is not so for some observables different from the large 
deviation function (e.g. for average profiles measured at 
the end of the large deviation event; see Ref. [22j]). In 
our case, the following choice turns out to be convenient 



Jc> 



y G [1, N — 1] (bulk exchange) 



N 1 (C9) 
y = 0, N (boundary baths) 



That is, we measure the energy current flowing through 
the bulk of the system. Using this current definition 
and eq. (|C8[) . we may write the modified normalized 



dynamics U' c , c = Y c l Uc y c ex P[^Jc y c)> which for y 6 
[l,iV-l] reads 



e A(e y -e' H ) 

U kc = YcW+T) 



(CIO) 



with A = \/(N - 1), while U' c , c = Y c 1 U C 'c for y 
0,N, see eq. (|C9|) . The exit rate is given by 



Y r , 



V 1 Ae„ „-Ae H+1 



N 



+ 1 + S A(iV + l)(e„ + e„ +1 )' (CU) 



In these paper we simulate a system of size N — 50, with 
Tl = 2 and Tr — 1, using M = 10 3 copies of the sys- 
tem and a maximum time of t = 10 Monte Carlo steps. 
For a given initial condition, we average the measured 
//(A) for different times once in the steady state, after 
a relaxation time of 2 x 10 3 Monte Carlo steps. In ad- 
dition, we average results over many independent initial 
conditions, in which local initial energies are randomly 
drawn according to the Gibbs distribution with temper- 
ature parameter T st [x — i/(N + 1)] corresponding to the 
linear, steady temperature profile. Fig 1 171 shows the con- 
vergence of /i(A) in time for a given value of A and many 
different initial conditions. Using the above method, we 
obtained an accurate measurement of the current large 
deviation function, see Fig. [3] in Section ITVl 



Appendix D: Time Reversibility and Statistics 
during a Large Fluctuation 

In this Appendix we use the time reversibility of 
the underlying stochastic dynamics to study the system 
statistics during a large deviation event and the symme- 
tries of the large deviation function and the associated 
optimal profiles, using the formalism described in Ap- 
pendix [C] In particular, we describe a relation between 
system statistics at the end of the large deviation event 
and for intermediate times. First, consider the probabil- 
ity P(Ct, Qt, t) that the system is in configuration Ct at 
time t with a total time-integrated current Q t . As in the 
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FIG. 17: (Color online) Time evolution of n(X) for A = —0.1 
and many different initial conditions. Here N = 50, M = 10 3 , 
and T L = 2, Tr = 1. 

previous appendix, we drop the dependence of this prob- 
ability on the initial state Co, which we assume lost for 
long enough times. This probability obeys the following 
master equation 

P(C t ,Q t , t)=J2 U Ct c'P(C, Q t - J Ct c>,t - 1) . (Dl) 
c 

which by iterating in time leads to 

t-i 

P(C u Q u t)= ]T U Ct c t - 1 ..U Cl c S(Qt-Y,Jc k+1 c k ), 

Ct-i-.Cx k=0 

(D2) 

and it is clear that P(Q t ,t) — J2c t P(Ct>Qt,t), see 
eq. (|C1|) in the previous appendix. Now, P™ d (Ct) = 
P(C t ,Qt,t)/P(Q t ,t;Co) is the probability of having a 



configuration Ct at the end of a large deviation event as- 
sociated to a current q = Qt/t. Defining H(C t ,X,t) = 
J2 Qt exp(XQ t )P(C t ,Qt, t) so that 

IL(C t ,X,t) = Yl U Ct c t _ 1 ...U Cl c , (D3) 

Ct— 

with Uc'C'W = Ucc ex P(^^cc): one can easily show 
that, for long times t, Px nd (C t ) = U(C u X,t)/Il(X,t) = 
Pq(y.(Ct), where q (^) is the current conjugated to pa- 
rameter A, and H(X,t) is defined in eq. (|C3I) . Using the 
spectral decomposition of Appendix [CJ it is simple to 
show that F g ond (C t ) cx {C t \A R (X)), so the right eigenvec- 
tor |A fl (A)) associated to the largest eigenvalue of matrix 
U(X) gives the probability of having any configuration 
at the end of the large deviation event. Noticing that, 
for the Monte Carlo algorithm described in the previous 
appendix, the fraction of clones or copies in state Ct is 
proportional to (C t \A H (X)) for long times, see eq. (|C7I) . 
we deduce that the the average profile among the set of 
clones yields the mean temperature profile at the end of 
the large deviation event, X^ nd (x). 

The initial and final time regimes during a large de- 
viation event show transient behavior which differs from 
the behavior in the bulk of the large deviation event, i.e. 
for intermediate times @. In particular, as we will show 
here, midtime and endtime statistics are different, though 
intimately related as a result of the time reversibility 
of the microscopic dynamics. Let P(C T , X,r,t) be the 
probability that the system was in configuration C T at 
time r when at time t the total integrated current is Qt- 
Timescales are such that 1 <C r <C t, so all times involved 
are long enough for the memory of the initial state Cq to 
be lost. We can write now 



P(C T) Q U T,t) = 



u, 



CtCt- 



• ••[/, 



C T + 1 C T ^ C tCt-1 



•U, 



CiCq 



>(Qt 



t-i 

E 

fc=0 



Jc k+1 c k 



(D4) 



where we do not sum over C T . Defining the 
moment-generating function of the above distribution, 
n(C T ,A, r,t) = Y,Q t exp(XQt)P{C T ,Q t ,T,t), we can 
again check that the probability weight of configuration 
C T at intermediate time r in a large deviation event of 
current q = Q t /t, Pf A {C T ) = P(C T ,Q U T,t)/P(Q t ,t), is 
also given by P{ nid (C T ) = IL(C T , A, r, t)/II(A, t) for long 
times such that 1 <C r -C t, with q — q a (X). In this 
long-time limit one thus finds 

Pr d (C T ) cx <A L (A)|CV)(CV|A*(A)) , (D5) 

in contrast to P™ d (C), which is proportional to 



(CIA^A)}, see above. Here |A fl (A)} and (A L (A)| arc the 
right and left eigenvectors associated to the largest eigen- 
value e A ( A ) of modified transition rate U(X), respectively. 
They are different because U is not symmetric. In order 
to compute the left eigenvector, notice that |A L (A)) is 
the right eigenvector of the transpose matrix U T (X) with 
eigenvalue e A ( x \ This right eigenvector of U T (X) can be 
in turn related to the corresponding right eigenvector of 
U(— A — £) by noticing that the local detailed balance 
condition holds for the KMP model, guaranteeing the 
time reversibility of microscopic dynamics. This condi- 
tion states that UccPeq(C) = UccPcq(C')e £ J cc j where 
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Peq(C) is an effective equilibrium weight which for the 
KMP model takes the value p cq {C) — exp(— Yl y =i Py e v) 
with C = {e y ,y = 1...N} and f3 y = T~ 1 + £§^ I . Local 
detailed balance then implies a symmetry between the 
forward modified dynamics for a current fluctuation and 
the time-reversed modified dynamics for the negative cur- 
rent fluctuation, i.e. Ucc — Pcq(C)U(—\ — £)p e<l (C), 
or in matrix form 

U T (X) = P-^C-A - £ )P eq , (D6) 

where P oq is a diagonal matrix with entries p eq (C). 
Eq. (|D6[) implies that all eigenvalues of U(X) and 
U{— A — £) are equal, and in particular the largest, so 
jLt(A) = fj,(— A — £) and this proves the Gallavotti-Cohen 
fluctuation relation. Moreover, if |A?(— A — £)) is a right 

eigenvector of U{— A — £), which can be expanded as 
|Aj*(-A-£)> = £ c (G|Af(-A-£))|C7), then 

c 

is the right eigenvector of t/ T (A) associated to the same 
eigenvalue. In this way, by plugging this into eq. (|D8|) 
we find 

Pf\C) (x (p e c ?)- 1 (C|A fl (-A-^))(C|A ii (A)) , 



where we assumed real components for the eigenvectors 
associated to the largest eigenvalue. Equivalently 

Pf ld (C)P cn A d JC) 
P A mid (C) = A 1 ; G q A ~ £V ; , (D8) 
Pc 

with A a normalization constant. This relation implies 
that configurations with a significant contribution to the 
average profile at intermediate times are those with an 
important probabilistic weight at the end of both the 
large deviation event and its time-reversed process. Sup- 
plementing the above relation with a local equilibrium 
hypothesis, one can obtain average temperature profiles 
at intermediate times in terms of profile statistics at the 
end of the large deviation event. 
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